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ABSTRACT 

Very- high- energy (VHE, E > 100 GeV) gamma rays emitted by extragalactic sources, such as 
blazars, initiate electromagnetic cascades in the intergalactic medium. The cascade photons arrive 
at the earth with angular and temporal distributions correlated with the extragalactic magnetic field 
(EGMF). We have developed a new semi-analytical model of the cascade properties which is more 
accurate than previous analytic approaches and faster than full Monte Carlo simulations. Within 
its range of applicability, our model can quickly generate cascade spectra for a variety of source 
emission models, EGMF strengths, and assumptions about the source livetime. In this Letter, we 
describe the properties of the model and demonstrate its utility by exploring the gamma-ray emission 
from the blazar RGB J0710+591. In particular, we predict, under various scenarios, the VHE and 
high-energy (HE, 100 MeV < E < 300 GeV) fluxes detectable with the VERITAS and Fermi Large 
Area Telescope (LAT) observatories. We then develop a systematic framework for comparing the 
predictions to published results, obtaining constraints on the EGMF strength. At a confidence level 
of 95%, we find the lower limit on the EGMF strength to be ~ 2 x 10 -16 Gauss if no limit is placed 
on the livetime of the source or ~ 3 x 10~ 18 Gauss if the source livetime is limited to the past ~ 3 
years during which Fermi observations have taken place. 

Subject headings: astroparticle physics — BL Lacertae objects: individual (RGB J0710+591) — cosmic 
background radiation — gamma rays: general — intergalactic medium — magnetic 
fields 



1. INTRODUCTION 

The extragalactic magnetic field (EGMF) is of great 
interest to the overall understanding of astrophysical 
magnetic fields and related processes. It could act 
as a see ding field for m agnetic fields in galaxies and 
clusters (jWidrowl I2002D . and its origin may be re- 
lated to inflation o r other periods in the ea rly his- 
tory of the universe (jGrasso fc Rubinstein! 1200 lh. Fara- 
day rotation measurements (iKronberg fc Perrvi Il982t 
lKronberd [l99l iBlasi et al.lll999ft and analysis of C OBE 
data anisotropy ( Barrow et al.lll997t iDurrer et al.ll2000h 
have put an upper bound on the EGMF strength 
at ~ 10~ 9 Gauss. On the other hand, gamma- 
ray-initiated electromagnetic cascades deflected by the 
EGMF in i ntergalactic space have a characteristic angu- 
lar sp read (|Aharonian et al.lll994f) and time delay (|Plagal 
Il995h . bot h of which provide a probe for lower EGMF 
strengths (|Neronov fc Semikozl 120071: lElviv et al.l I2009D . 
Present and next generation gamma-ray telescopes have 
the possibility to measure the EGMF strength by observ- 
ing the angular and temporal distributions of cascade 
photons from extragalactic gamma-ray sources such as 
blazars (|Neronov fe SemikozfeoOl iDolag et al.l[2009h . 

With current very-high-energy (VHE, E > 100 GeV) 
and high-energy (HE, 100 MeV < E < 300 GeV) 
gamma-ray data on VHE-selected blazars, a lower limit 
on the EGMF strength can be placed by requiring that 
the cascade flux of VHE emission not ex ceed the mea- 
sured flux or upper lim i t in t he HE band ( Muras e~et al.l 
120081 INeronov fc Vovkl I2010D . Analytic cascade mod- 
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els assuming a simple relationship between the cascade 
flux and EGMF strength have put the lower limit at 
10~ 16 to 10~ 15 Gauss when the source livetime is un- 
limite d (|Tavecchio et al.l l201Ct 120111: INeronov fc Vovkl 
120101 ), similar to the results of Monte C arlo simula- 
tions (|Dolag et al.ll2011l : lTavlor et al.ll2011[ L If the cas- 
cade time delay is limited to the ~ 3 years of simultane- 
ous HE and VHE observations, the lower limit becomes 
10 -19 t o 10~ 18 Gauss accor ding to the simple cascade 
models (|Dermer et all 120111) or 10~ 18 to 10~ 1 7 Gauss 
according to the simulations (jTavlor et al.ll20li[) . 

In this Letter, we present a new semi-analytic model 
of the electromagnetic cascade. In contrast to pre- 
vious analytic cas c ade models (iDermer "eTall I20T1I: 
INeronov fc Semikozl 120091 : iTavecchio et al.l l2011h . our 
cascade model considers the full track of the primary 
photon without the assumption of interacting exclusively 
at the mean free path. This simultaneously accounts for 
both angular and temporal constraints in a natural way. 
In addition, we model the radiation backgrounds and 
source emission in greater detail. As a co mplementary 
approach to full Monte Carlo simulations (jDolag et al.l 
12011 iTavlor et al.1 MHi lArlen et al.ll20ll . our model 
serves as a tool for clarifying the cascade picture, rapidly 
searching through the parameter space, and interpreting 
simulation results. We also develop a systematic frame- 
work for applying our cascade model's predictions to de- 
rive lower limits on the EGMF strength at specific confi- 
dence levels. This framework is applicable to the results 
of Monte Carlo simulations as well. 

2. MATHEMATICAL MODEL OF CASCADES 

Consider a blazar emitting gamma rays at distance L. 
Gamma rays emitted at an angle 9 S relative to the line of 
sight may produce an e ± pair via absorption on the pho- 
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Fig. 1 . — Geometry used t o calculate the gamma-ray cascade, 
following IDermer "eTatl j201ll V 



ton background ([Gould fc Schrederlll967ft at a distance 
L' from the source, not necessarily equal to the mean 
free path. After being deflected by the EGMF through 
the angle 9d, the pairs could scatter background photons 
to gamma-ray energies, redirecting them toward the ob- 
server. These secondary gamma rays would arrive at an 
incidence angle 9 C (see Fig. [TJ). In this picture, the influ- 
ence of the EGMF enters solely through 9d- The angles 
9 C and 9 S are uniquely specified in terms of Q4, L, and 
L', provided 9 C < n/2. 

Because the energy density of CMB photons far ex- 
ceeds that of the extragalactic background light (EBL), 
we assume that inverse Compton scattering proceeds in 
the Thomson regime via e^- interactions with the CMB 
exclusively. An electron with Lorentz factor 7 e will on 
average scatter secondary photons to energy 47 2 eo/3, 
wher e ep ~ 0.64 meV is the ave rage CMB photon en- 
ergy (jBlumenthal fc Gouldlll970() . The energy loss rate 
of the electron is 



dj e m e c 2 
dt 



-g£ocncMBCTT7 e 



(1) 



where c is the speed of light, uqmb ~ 411 cm -3 is the 
CMB photon density, and ut ~ 6.65 x 10 -25 cm 2 is the 
Thomson cross section. For an EGMF B perpendicu- 
lar to the electron momentum p e , the deflection rate in 
terms of the Larmor radius 77 — 7 e m e c 2 /ei? is 



d9d 
dt 



c _ eB 
n j e m e c 



(2) 



Thus, when an electron's Lorentz factor has changed 
from 7 e o to j e , the electron will have been deflected by 

N 3 eB ( 1 1 \ 

0do(7eO,7e) = o — • 3 

8eo«CMBCTT \li 7e0/ 

However, if the angle 9f between B and p e is other than 
7r/2, Eq. El should be generalized to 

0d (7e0 > 7e> Of) = arccos (sin 2 9 f cos 9 d0 + cos 2 6 } ) . (4) 

We use the full CMB black-body spectrum to deter- 
mine the observed cascade spectrum. The number of 
secondary photons between energies E and E + dE pro- 
duced by an electron that changes Lorentz factor from 
7e + c?7 e to 7 e is the number density of CMB photons be- 
tween 3S/(47 2 ) and 3(E+dE)/(4j 2 ) times the Thomson 
cross section and the distance travelled by the electron: 



dN(E,j e )=cdta T 
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dE 



8 7e 4 /j3 c 3 7 2( e 3E/4( 7 2fcT) _ 1) 

81irE 2 m e d'y e dE 
32/i 3 c7 e 8 e n C MB e^K^T) _ 1 



(5) 



where h is the Planck constant, k the Boltzmann con- 
stant, and T = 2.73 K the CMB temperature. 

In terms of the mean free path A(e) = L/r{e) inferred 
from the optical depth r(e) of the EBL, the probabil- 
ity for a photon of energy e to be absorbed between 
V and L' + dL' is e~ L '> dL' / 'A(e) . Approximating 
both particles in the resulting pair to have initial energy 
TO e 7eoc 2 = e/2, we calculate the differential flux of ob- 
served secondary photons by integrating over L' , e, and 
7 e , and averaging over 9 / : 

d r y e 81irE 2 m e 



dN(E) 
dE 



d9 f g(9 f ) J dej dL' 



16h 3 c^e n C MB (e 3£ /( 4 ^ fcT ) - l) 



A(e) 



-f(e,e s ) 



cxp i-^jL 2 + L' 2 -2LL' 



d/KE)) (6) 



Here, g(9f) is the probability distribution of 9f, equal 
to sm Of for an EGMF uniformly distributed in direc- 
tion, and f(e,9 s ) is the intrinsic flux of the source, with 
9 S = d - 9 C from Fig. [Q The 1/2tt factor from the e ± 
being deflected into the surface of a cone with opening 
angle 9d cancels the 2-7T enhancement from a similar ef- 
fect at the source. We take the integral over Of from 
to 7r/2. The physical lower bound for the integration 
over the primary energy e is e m i n = 2 r y e m e c 2 . Nearly all 
of the photons above 200 TeV will be absorbed within 1 
Mpc of the source, and the e ± pairs will be isotropized 
by the strong field of the surrounding galaxy, resulting 
in negligible cascade contribution. We therefore adopt 
an upper limit of e max = 200 TeV, suggesting an upper 
limit of e max /2m e c 2 on the 7 e integration. As a practi- 
cal matter, we enforce a lower limit of j e = 10 on the 
7 e integration, motivated by the CMB density becoming 
negligible at energies above 3 meV and our disinterest in 
the cascade spectrum below 100 MeV. 

Observational effects enter Eq. [6] through limits on the 
L' integration. As seen in Fig. [TJ we can express C as 



iircsiii I — s'mOd 



(7) 



so that a cut on 9 C translates directly into a cut on L'. 
Similarly, the time delay AT of cascade photons may be 
written as 



cAT = L' + ^L 2 + L' 2 - 2LL' cos{9 d - 9 C ) - L, (8) 

exchanging a limit on the source livetime AT for another 
constraint on L' . We adopt the intersection of the L 1 cuts 
from Eqs.[7]and[S]in evaluating the cascade flux via Eq.[5J 
We now briefly examine several assumptions made in 
the construction of this cascade model, (i) Exact energy 
distributions of pair production products are approxi- 
mated as each having half the energy of the primary 
photon, as the cascade spectrum o nly weakly depends 
on th e spectral distribution of pairs (jCoppi &: Aharonianl 
1997). (ii) The Thomson limit assumption, appropriate 
when \JTqE <C m e c 2 E <C 400 TeV, certainly produces 
spectra valid for the range of existing TeV gamma-ray 
instruments, (iii) We assume the EGMF to be coher- 
ent over the cooling length of the electrons, at most a 
few Mpc for the most energetic electrons, making this 
assumption valid for coherence lengths > 1 Mpc. 
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(iv) Only secondary cascade photons are considered. For 
observed photons above 100 MeV, the primary energy 
must be e > 2m e c 2 y / 3 x (100 MeV)/4e « 400 GeV. 
The requirement for third-generation photons is thus e > 

2m e cV 3 x ( 40 ° GeV)/4e w 25 TeV. While a primary 
photon above this energy does produce higher-generation 
cascade photons, the power going into the higher gen- 
erations is small for conventional blazar emission mod- 
els, leadin g to negligible contrib ution from higher-order 
cascades (|Tavecchio et al.l 1201 lh . (v) Cosmological ef- 
fects enter the cascade model solely in the calculation 
of A(e). Cosmic expansion, energy redshift, EBL evolu- 
tion, and other cosmological effects are ignored, limiting 
the application of the cascade model to nearby sources 
(z < 0.2), considering the (1 + z) 4 radiation density evo- 
lution, (vi) We assume axial symmetry in the intrin- 
sic emission /(e, 9 S ), following the current understanding 
of blazar emission, that is, boosted isotropic emission 
from hot blobs with the jet dire ction pointing toward 
the earth (jUrrv fc Padovariilll995l ). (vii) The blazar in- 
trinsic emission should be steady over the time interval 
AT. 

To get the predicted cascade flux from Eq. [6l 
we require a source emission model f(e,9 s ). Moti- 
vated b y models of blazars as relativistically beamed 
sources (|Urrv fc PadovanilH995h we model the intrinsic 
emission as boosted isotropic emission, 



- intrinsic Spectrur 

- B = 10" 17 Gauss 

- B = 10" 16 Gauss 

- B = 10" 1S Gauss 



,r a - l e- c 



e/E 



f(e, 9 S ) = f (1-/8 cos ( 

+ fo(l + Pcos9 s )- a - L e- a e^*°, (9) 

that is, a power law of index a boosted by the Lorentz 
factor T = 1/ ^/l — f3 2 with an exponential cutoff en- 
ergy Eq. The second term in Eq. [9] models the blazar 
counter jet, which we find does not significantly impact 
our result. To get a conservative esti mate for the cascade 
flux, w e employ the EBL model from iFranceschini et all 
(2008) which is relatively transparent for VHE gamma 
rays. Having no prior assumptions on the EGMF struc- 
ture, we take g(0f) = sm9f to calculate the cascade flux 
in Eq. [6l 

3. APPLICATION OF THE MODEL TO BLAZAR DATA 

We demonstrate the cascade model by investigating 
RGB J0710+591, a high-frequency-peaked BL Lacer- 
tae (HBL) object located at a redshift of z = 0.125, 
for which simultaneous HE and VHE data are avai l- 
able with no variability observed (jAcciari 
In the VHE regime, we take the spectrum reported 
by lAcciari et al.l ([2010), while to get the HE spectrum 
we analyze publicly available Fermi Large Area Tele- 
scope (LAT) data taken from a ~ 3-year period between 
2008 August and 2011 January, using the Fermi Science 
Tools v9rl8p6 and the P6_V3_DIFFUSE instrument re- 
sponse functions (IRFs), with models gll_iem_v02 for the 
galactic diffuse and isotropic_iem_v02 for the isotropic 
background 3 . Accounting for nearby point sources, we 
perform an unbinned likelihood analysis to get the spec- 
trum between 100 MeV and 300 GeV, finding a best-fit 
index of 1.62 ± 0.11. Next, we bin the data into five en- 
ergy bins, demanding that each bin have a test statistic 
greater than 10 and a maximum relative uncertainty on 

3 http:/ /fermi. gsfc.nasa.gov/ssc/ 
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Fig. 2. — Measured spectrum of RGB J0710+591 (points) and 
predictions for different situations (curves). Squares: VERITAS 
data points. Shaded area: Fermi 68% confidence band. Circles: 
Fermi data points. Panel a: Total fluxes for a source with un- 
limited livetime. Panel b: Total fluxes for a source with 3-year 
livetime. 



the flux of 50%, as suggested by lAbdo et all (|2010[ ). The 
combined spectra appear in Fig. [21 demonstrating that 
the Fermi spectral points are consistent with the 68% 
confidence band from the unbinned likelihood analysis. 

To compare with the measured spectra, we compute 
the total flux as the sum of the cascade flux calculated 
by Eq. [5] and the direct flux /(e, 0) exp (— r(e)), and then 
normalize it to observational data. For example, the re- 
sults for a = 1.5, Eq = 25 TeV, T = 10, and a range 
of different EGMF strengths appear in Fig. [Ha) for a 
blazar with unlimited livetime, and in Fig. [^b) for a 
blazar active for only 3 years. In both cases we require 
9 C to be smaller than the 68 % containment radius of the 
point spread function (PSF) (|Atwood et al.ll2009D in the 
Ferm i range, or smaller than the 9 2 cut ( Acciari et al.l 
120101) in the VHE range. Setting the simple require- 
ment that the flux be lower than the measured spec- 
trum, we observe from Fig. [5] that the lower limit on 
the EGMF strength should be between 10" 16 and 10~ 15 
Gauss for the infinite-time case, and between 10 -18 and 
10~ 17 Ga uss for the 3-year c ase, consistent with results 
shown bv lTavlor et all (|20Tll) for RGB J0710+591. 

To be more rigorous, we fit the predicted flux to the 
data and use the resulting \ 2 values to constrain the 
EGMF strength. Fig.[3fa) shows a sample fit for a = 1.5, 
Eq = 25 TeV, V = 10, and B = 3 x 10~ 16 Gauss. 
Because we cannot exclude the existence of additional 
components contribu ting to the blazar emission (e.g., 
iBottcher et al.l [20081) and modifying the intrinsic spec- 
trum of Eq. [JU we also fit a broken power law with an 
index abrcak below 80 GeV, as shown for instance in 
Fig. EJb). We find that the additional free parameter 
does not greatly affect our constraints on the EGMF. 

For each EGMF strength B, we find the best-fit model 
of the combined cascade and intrinsic emission for a wide 
range of cutoff energies Eq G (0.1 TeV, 100 TeV), allow- 
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Fig. 3. — Fitting RGB J0710+591 spectrum with cascade model 
predictions assuming unlimited source livetime. Squares: VERI- 
TAS data points. Shaded area: Fermi 68% confidence band. Cir- 
cles: Fermi data points. Panel a: Fitting the sum of cascade flux 
and direct flux. Panel b: As Panel a, except the direct flux power- 
law index changes below 80 GeV (see text). 



12 



10 





I I I 

— 1 year 
10 years 

- - - 100 years 
■"Vtf * — 10 3 years 

\ '"X."'^"'""'' y ears 

^ \\ \ 10 5 years — 

\ \i \ ■ ■ ■ 1 6 years 
'* ^ \"s^» — Unlimited 

\ \ \ ► 

\ \ \ : 
\ \ \ i 






• 90% constraint \ 

_ ■ 95% constraint \^/ 


( x X i ssf 




1 1 1 M 



10" 1 



10" 1B 10~ 17 10"' 

EGMF Strength [Gauss] 



10" 1 



10' 1 



Fig. 4. — Best-fit x 2 versus EGMF strength B for various live- 
times for RGB J0710+591. Lower limits on B are indicated at 90% 
and 95% confidence levels for the 1-year and unlimited cases but 
omitted at intermediate times for clarity. 



ing a and abreak to vary between 2.5 and the physically 
motivated constraint of 1.5 (e.g., iMalkov fc O'C Drurvl 
1200 It lAharonian et al.l I2006D and fixing T to a typical 
value of 10. We plot the x 2 value of the best-fit model as 
a function of B in Fig. 2] for blazar livetimes from 1 year 
to 10 6 years, and for the unlimited case. At low B, the 
X 2 values converge because the cascade arrives promptly 
and at small angles. As B increases, the arrival angles 
and times of the cascade begin to spread out, diminish- 
ing the observed emission and providing a better fit to 
the observed data. 

The convergence of the curves in Fig. 0] to the infinite- 
livetime curve is easily understood. Combining Eq. [7] 
with Eq. [8] and assuming small angles, we can translate 
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Fig. 5. — Excluded field strength as a function of blazar livetime 
for two different confidence levels for RGB J0710+591. The solid 
black lines of constant B 2 JT agree with the curves for small times 
but begin to diverge for livetimes between 10 4 and 10 5 years, as 
expected. At longer times, the instrument PSF limits the possible 
delay time of the photons, causing the constraint to level off. 



a cut on the instrument angle C into a time cut AT f 

(L-L')L q2 
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For example, source photons at 1 TeV, which have a 
mean free path of L' ss 400 Mpc, will produce cascade 
photons of energy E s» 0.8 GeV, for which an angular 
cut of 8 C ~ 1° is appropriate for the Fermi LAT. At 
the distance of RGB J0710+591 (L w 500 Mpc), this 
translates into a time cut of AT C w 6 x 10 4 years. This 
is close the the livetime of ~ 10 5 years at which the 
X 2 curves in Fig. Q] begin to converge to the unlimited- 
livetime case, becoming nearly indistinguishable at <~ 10 6 
years. If the blazar livetime is smaller, the livetime cut 
outweighs the angle cut, and the position of the \ 2 curve 
depends on the blazar livetime. For longer livetimes, 
the angle cut becomes more constraining and the curves 
converge to the unlimited-livetime case. 

We constrain the EGMF strength for a given livetime 
by finding the point at which \ 2 exceeds its minimum 
value by A% 2 for each curve in Fig. [H Two sample con- 

3Spondin 
3sf l2006|) 

in the figure. The limits derived from these confidence 
levels are shown in Fig. [5] as a function of the blazar 
livetime. For livetimes below ~ 10 4 years, the limit on 
the EGMF strength scales with the blazar livetime AT 
as B ~ v 7 AT, as expected from combining Eq. [2] and 
Eq. 1101 This relation breaks down when the time limit 
enforced by Eq. [10] is smaller than the blazar livetime. 
Our limit at 95% confidence is B > 2 x 10 -16 Gauss if 
the blazar's livetime is infinite and B > 3 x 10~ 18 Gauss 
if the livetime i s 3 ye ars, in agreement with the results 
of I Taylor et al.l (|201 lh - In all cases, we reject the hypoth- 
esis of zero EGMF at greater than 2<r confidence. 

4. DISCUSSION AND CONCLUSION 

The cascade model we present builds upon previous 
models to give a more complete semi-analytic treatment, 
accounting for the photon trajectories, background spec- 
tra, and intrinsic emission in sufficient detail to produce 
realistic cascade predictions for different blazar livetimes. 



fidence levels (90% and 95%), correspondi ng to Ax 2 val- 
ues of 2.72 and 3.84 (see e.g., Uamesil2006D are indicated 
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For example, the spectra shown in Fig. [2] converge to 
the zero field case at high energies and become sup- 
pressed by the field approximately as 1 /B at low energies 
(~ 1 GeV) , in agreement with the spectra p roduced by a 
full M onte Carlo simulation as pre sented inlTavlor et all 
(1201 lh . Previous ca scade models (jTavecchio et alJ[2010t 
iDermer et alj[201ll) overestimate the suppression of the 
cascade flux in the EGMF compared to our result. 

Combined with a systematic framework for interpret- 
ing the simultaneous fit of HE and VHE spectra to pre- 
dictions, our model derives robust lower limits on the 
EGMF strength. The application of the model to RGB 
J0710+591 data yields lower li mits at 95% confid ence 
that agree with the results of I Taylor et all (|2011|) for 
the same source, further indicating that the physical pic- 
ture of the cascade presented by the model is sufficiently 
complete to produce accurate results within the limits 
discussed in Sec. [2] The scaling of the lower limits with 
source livetime, depicted in Fig.0 demonstrates that our 
analysis framework is also physically meaningful. The 
limits are conservative because we choose a low EBL 
model, as well as for a few other reasons that we now 
discuss. 

We assume that the EGMF is configured as coherent 
domains of the same field strength and random field ori- 
entation, neglecting domain crossing by the pairs. In 
the energy range of interest, this restricts the domain 
size As to be As > 1 Mpc. If A# < 1 Mpc instead, do- 
main crossing will result in smaller deflection ang les for 
pairs (jlchiki et alj|2008t iNeronov fc Semikozll2009b and a 
higher bound on the field strength than we report. 

In testing various models of the intrinsic emission, we 
only consider cutoff energies below 100 TeV and spec- 

4 http: // fermi . gsfc. nasa.gov / ssc / data / analysis/L AT_caveats . html 



tral indices softer than 1.5. While it is possible to 
achieve intrinsic emission that is even har der by invok- 
ing s o me unconventional mechanism (e.g.. IStecker et all 
120071; iBottcher et "all 120081: lAharonian et all |2008| ). the 
HE cascade flux could only be higher in that case, leading 
to a more stringent lower limit. The spectral index could 
also be furthe r constrained using mu lti- wavelength infor- 
mation (e.g., iTavecchio et al.l [201 lh . which depends on 
detailed blazar modeling and hence is beyond the scope 
of this Letter. 

The P6.V3_DIFFUSE IRFs used in the cascade model 
integ ration may underestimate the PSF at large ener- 
gies (|Ando fc Kusenkol I2010t INeronov etaLl l2011h . but 
Fig. [2] shows that the most constraining part of the HE 
spectrum for RGB J0710+591 is the low-energy region 
where this effect is small , and, as discussed in Sec. O it 
can only affect the cases with livetimes larger than 10 4 
years. Furthermore, we expect the IRFs we use to un- 
derestimate the cascade flux, so our lower limit is still 
conservative. In the future, there will be more realistic 
IRFs for the Fermi LAT, as well as more blazars with 
simultaneous HE- VHE baseline data available, and the 
lower limit could be further improved using the frame- 
work we present. 
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